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STATEMENT  OF  PROBLEM  STUDIED 


Background 

The  Department  of  Defense  (DOD)  has  thousands  of  sites  with  subsurface  contamination. 
The  Army  is  charged  with  the  responsibility  of  coordinating  DOD’s  technical  efforts  toward 
cleanup  of  these  sites.  It  is  well  understood  that  efficient,  accurate  computer  models  of 
subsurface  flow  and  transport  are  vital  tools  in  site  characterization  and  in  assessment  of 
remediation  strategies.  Also  important,  but  perhaps  less  widely  realized,  is  the  role  that  such 
models  can  be  expected  to  play  in  fundamental  understanding  of  field-scale  phenomena,  by 
numerically  upscaling  smaller-scale  experimental  results  through  a  hierarchy  of  scales  of 
heterogeneity. 

The  objective  of  this  research  was  to  design  and  implement  new  numerical  techniques  that 
will  be  particularly  efficient  and  accurate  for  simulation  of  groundwater  flow  and  transport. 
The  prototypical  problem  considered  was  two-  and  three-dimensional  solute  transport  in  the 
saturated  zone,  to  be  followed  in  future  research  by  extensions  to  vadose-zone  transport  and 
multiphase  flow.  Simultaneously,  model  development  efforts  with  more  conventional  methods 
were  carried  forward  by  the  staff  at  the  Waterways  Experiment  Station  (WES).  The  research 
of  this  project  offers  the  potential  of  substantially  greater  efficiency  and  accuracy,  but  with 
less  certainty  of  success  owing  to  the  newness  of  the  approaches.  The  products  are  theoretical 
understanding  of  and  practical  experience  with  faster,  better  methods  for  flow  and  transport, 
and  more  tangibly,  code  modules  (still  under  development)  fitting  into  the  structure  of  the 
WES  model  that  will  offer  the  new  techniques  as  choices  to  the  user. 

Field  groundwater  modeling  problems  are  physically  three-dimensional  and  usually  ex¬ 
hibit  significant  heterogeneity  in  hydraulic  conductivity,  porosity,  and  other  properties.  Nev¬ 
ertheless,  models  often  assume  that  a  lower-dimensional  or  more  homogeneous  representation 
is  adequate.  This  is  frequently  the  result  of  practical  constraints,  namely  the  inability  of 
available  software  to  solve  the  real  problem  accurately  within  available  computing  resources. 
One  would  like  to  come  closer  to  a  state  of  affairs  in  which  decisions  to  simplify  a  model 
could  be  based  on  scientific  validity  rather  than  practical  necessity. 

Transport 

Difficulties  with  Eulerian  methods.  Serious  troubles  occur  in  approximations  of  transport 
equations  when  advection  dominates  diffusion  and  dispersion.  Standard  Eulerian  numerical 
methods,  such  as  centered  finite  differences  or  Galerkin  finite  elements,  perform  well  when 
the  grid  size  is  sufficiently  small.  In  many  practical  situations,  this  is  not  feasible,  and  such 
methods  produce  nonphysical  oscillations  in  concentration.  The  usual  remedy  is  to  remove 
them  by  employing  an  upstream  finite  difference  or  finite  element  method.  Resulting  con¬ 
centration  profiles  do  not  oscillate,  but  they  are  seriously  dispersed  relative  to  their  physical 
counterparts.  The  dilemma  is  a  choice  between  spurious  oscillations  and  numerical  disper¬ 
sion.  In  practical  nonlinear  problems,  oscillations  will  cause  serious  distortions  and  cannot 
be  allowed  if  accuracy  is  to  be  maintained.  The  same  is  true  of  numerical  dispersion.  Addi¬ 
tionally,  one  of  the  principal  objectives  in  fundamental  understanding  of  field-scale  transport 
is  to  properly  upscale  dispersive  phenomena  to  account  for  a  hierarchy  of  heterogeneities; 
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a  model  is  clearly  useless  in  this  regard  if  it  overwhelms  physical  dispersion  with  numerical 
dispersion.  The  dilemma  is  not  the  result  of  the  physical  process;  the  grid  needed  to  avoid 
oscillations  is  at  a  finer  scale  than  that  of  the  traveling  advection  fronts,  and  other  methods 
can  be  more  efficient. 

Advection  domination  also  leads  to  large  time-truncation  errors  for  Eulerian  methods. 
When  a  steep  front  passes  by,  as  typically  happens  when  advection  dominates,  errors  will 
be  large  unless  time  steps  are  very  small.  A  more  appropriate  method  would  reduce  these 
errors  by  time  stepping  in  a  Lagrangian  manner  that  would  follow  the  flow,  permitting  larger 
but  still  accurate  time  steps.  Thus,  the  methods  on  which  almost  all  existing  public-domain 
transport  codes  are  based  are  ill-suited  for  advection-dominated  transport,  both  spatially 
and  temporally. 

Lagrangian  methods.  Pure  advection  problems,  in  which  diffusion  and  dispersion  are  absent, 
can  be  solved  efficiently  and  accurately  by  Lagrangian  schemes  based  on  the  method  of  char¬ 
acteristics.  Such  methods  normally  build  on  propagation  of  discontinuities  (sharp  fronts), 
explicitly  tracking  front  positions.  Being  explicit,  these  methods  usually  require  that  the 
Courant  stability  condition  be  satisfied,  and  they  are  difficult  to  adapt  to  problems  with  sig¬ 
nificant  dispersion.  For  subsurface  contaminant  flow  and  transport,  fronts  can  be  steep  but 
will  not  be  discontinuous  because  of  dispersion,  and  advection-  and  dispersion-dominated 
regimes  will  both  be  important,  often  in  different  portions  of  the  same  computation.  The 
research  of  this  project  combined  Eulerian  and  Lagrangian  concepts  in  order  to  obtain  the 
advantages  of  both,  as  described  in  the  results  section. 

Flow 

Difficulties  with  anisotropy  and  heterogeneity.  The  principal  challenge  with  flow  equations 
is  to  obtain  accurate  fluid  velocity  fields  (or,  equivalently,  accurate  streamlines)  when  the 
subsurface  formation  is  highly  anisotropic  and/or  heterogeneous  in  its  hydraulic  conductivity, 
and/or  geometrically  irregular  in  its  geology.  This  situation  is  typical  in  field  cases,  where 
a  hierarchy  of  length  scales  of  heterogeneity  can  be  expected.  Since  desirable  methods  for 
transport  should  have  a  Lagrangian  component,  accuracy  in  transport  calculations  depends 
on  accurate  flow  velocities  as  input  to  the  transport  methods. 

At  interfaces  between  geological  layers  or  other  features,  large  spatial  discontinuities 
in  the  hydraulic  conductivity  are  routine.  Across  such  an  interface,  the  pressure  gradient 
is  discontinuous,  but  the  velocity  or  flux  is  comparatively  well-behaved,  having  continuous 
normal  component.  Thus,  procedures  that  solve  for  the  flux  directly,  and  whose  errors  depend 
on  appropriate  smoothness  of  the  flux  rather  than  the  pressure,  would  be  advantageous. 
Irregular  (non-rectangular)  geology  adds  another  layer  of  complexity  that  typical  models  do 
not  treat  rigorously. 
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SUMMARY  OF  THE  MOST  IMPORTANT  RESULTS 


Formulation  and  implementation  of  transport  codes 

ELLAM.  For  transport  equations,  the  project  studied  Eulerian-Lagrangian  localized  adjoint 
methods  (ELLAM).  These  methods,  which  are  based  on  space-time  finite  elements  whose 
temporal  faces  follow  space-time  streamlines,  generalize  earlier  Eulerian-Lagrangian  meth¬ 
ods  and  overcome  their  disadvantages.  A  major  one  is  that  such  procedures  are  not  mass- 
conservative,  and  though  the  errors  are  usually  small  they  sometimes  grow  to  troublesome 
magnitudes  greater  than  1%.  In  contaminant  transport  calculations,  where  small  concen¬ 
trations  related  to  regulatory  levels  are  important,  such  mass  losses  cast  doubt  on  otherwise 
useful  results.  A  further  vexing  difficulty  is  the  treatment  of  backtracked  streamlines  that 
cross  an  inflow  boundary  under  certain  types  of  boundary  conditions;  it  is  not  at  all  clear 
how  to  even  formulate  earlier  Eulerian-Lagrangian  methods  under  some  such  circumstances. 

The  ELLAM  formulation  was  suggested  by  Celia,  Russell,  Herrera,  and  Ewing  in  1990. 
ELLAM  multiplies  a  conservation  equation  by  a  space-time  finite-element  test  function  and 
integrates  over  the  space-time  domain.  The  key  is  to  choose  this  test  function  to  be  con¬ 
stant,  or  nearly  so,  on  space-time  streamlines,  and  to  define  space-time  finite  elements  whose 
temporal  faces  follow  streamlines.  The  integral  formulation  is  naturally  mass-conservative 
and  adapts  to  any  boundary  condition  in  a  systematic  way.  Prior  to  this  project,  ELLAM 
had  been  implemented  and  tested  for  model  transport  equations  in  one  and  two  dimensions. 
On  the  theoretical  side,  for  relatively  simple  one-  and  two-dimensional  cases,  Reference  1 
contained  convergence  proofs  that  verified  that  ELLAM  will  preserve  the  efficiency  and  accu¬ 
racy  advantages  of  earlier  Eulerian-Lagrangian  methods.  These  proofs  account  for  all  types 
of  boundary  conditions  and  for  temporal  discretization  of  the  outflow  boundary. 

Code  development.  For  solute  transport  with  variable  velocity  field  in  two  dimensions,  Ref¬ 
erence  8  documented  an  implementation  in  the  USGS  code  VS2D.  This  version  is  a  control- 
volume  formulation,  with  piecewise-constant  test  functions,  that  has  local  conservation  prop¬ 
erties  and  is  well-suited  for  incorporation  in  existing  finite  difference  codes  (such  as  VS2D). 
It  has  performed  very  well  in  heterogeneous  linear  and  homogeneous  quarter  five-spot  tracer 
flows. 

The  accuracy  of  the  method  depends  on  the  accuracy  of  the  Lagrangian  tracking  algo¬ 
rithm  that  defines  the  space-time  streamlines.  The  tracking  algorithm  developed  for  this 
code  (References  12,  19,  26)  is  exact  for  velocity  fields  that  have  the  lowest-order  Raviart- 
Thomas  form  (^-component  is  continuous  piecewise-linear  in  x  and  t,  piecewise-constant  in 
y;  analogous  for  y-component).  The  flow-equation  methods  discussed  below  deliver  accurate 
velocities  of  precisely  this  form.  This  tracking  scheme  has  proved  to  be  particularly  impor¬ 
tant  in  velocity  fields  in  which  a  particle  can  reverse  direction,  where  other  tracking  methods 
are  prone  to  significant  errors. 

Source  terms  were  later  added  to  this  code  (References  14,  21).  Then  the  formulation  was 
extended  to  three  space  dimensions  and  implemented  in  the  USGS  method-of-characteristics 
package  MOC3D  (References  15,  22,  24,  29).  The  greatest  difficulties  lay  in  the  integration 
of  solute  concentrations  over  potentially  irregular  regions  of  three-dimensional  space.  The 
ELLAM  transport  code  was  linked  to  the  USGS  MODFLOW  flow  code,  and  tested  on 
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a  suite  of  sample  problems  developed  for  MOC3D.  The  code  obtains  accurate  results  with 
long  time  steps.  The  velocity  information  provided  by  the  flow  code  developed  in  this  project 
is  of  the  same  form  as  that  provided  by  MODFLOW,  so  the  link  to  the  new  flow  code  will 
be  a  straightforward  modification  of  the  link  to  MODFLOW.  The  input  data  structure  is 
compatible  with  MODFLOW  inputs. 

Computational  results.  In  some  cases  the  improvements  over  other  methods  were  spectacular; 
for  one  sample  problem,  ELLAM  obtained  in  seven  time  steps  a  solution  similar  to  those  that 
a  standard  finite-element  code  and  earlier  MOC3D  codes  needed  hundreds  or  thousands  of 
steps  to  compute.  To  summarize  the  results  briefly,  ELLAM  can  treat  advection-dominated 
transport  with  coarser  grids  and  longer  time  steps  than  other  methods,  without  suffering 
from  nonphysical  oscillations  in  the  results.  Additionally,  ELLAM  results  are  relatively 
insensitive  to  the  orientation  of  the  grid  with  respect  to  the  direction  of  the  flow  velocity. 
This  means  that  one  can  obtain  accurate  simulations  on  simple  rectangular  grids  that  are 
much  easier  to  handle  than  irregularly-shaped  grids. 

Theory.  Advances  were  also  made  in  the  theoretical  convergence  analysis  of  the  control- 
volume  ELLAM.  It  can  be  viewed  as  a  Lagrangian  finite-volume  element  method,  and  the 
groundwork  for  analysis  of  time-dependent  Eulerian  methods  of  this  type  was  laid  in  Ref¬ 
erence  5.  Reference  25,  in  preparation,  is  the  first  extension  of  these  ideas  to  the  Eulerian- 
Lagrangian  framework. 

Formulation  and  implementation  of  flow  codes 

Mixed  finite  element  methods.  The  procedures  investigated  for  flow  equations  in  this  project 
were  mixed  finite  element  methods  (MFEM).  The  ELLAM  transport  algorithm  depends  for 
its  accuracy  on  accurate  streamlines.  As  shown  by  several  investigators,  the  MFEM  produces 
better  streamlines  in  groundwater  applications  than  standard  approaches,  particularly  when 
hydraulic  conductivity  is  heterogeneous. 

Raviart-Thomas  method.  The  first  MFEM  considered  in  this  work  was  introduced  by  Raviart 
and  Thomas  in  1976.  The  fundamental  idea  is  to  decompose  the  flow  equation  into  Darcy’s 
law  and  conservation  of  mass,  multiply  each  equation  by  its  own  test  function,  and  solve 
the  resulting  system  for  the  velocity  (flux)  and  the  pressure.  To  discretize  the  system, 
elements  (rectangles,  triangles,  or  parallelograms)  are  chosen  and  shape  and  test  functions 
are  specified.  This  is  done  in  such  a  way  that  the  velocity  vector  field  is  determined  by  its 
normal  fluxes,  which  are  continuous  across  edges  or  faces. 

The  advantages  of  this  approach  are  considerable.  The  velocity  field  conserves  mass  lo¬ 
cally,  so  that  it  can  be  used  in  the  ELLAM  transport  algorithm,  with  the  exact  tracking  used 
there,  without  creating  or  destroying  mass.  Another  major  advantage  is  that  the  accuracy 
of  the  approximate  velocity  can  be  shown  theoretically  to  depend  on  the  smoothness  of  the 
velocity,  not  that  of  the  pressure.  As  discussed  previously,  in  heterogeneous  media  this  favors 
MFEM  over  standard  methods  and  justifies  theoretically  the  computational  observations  of 
many  investigators,  including  some  in  which  flow  was  coupled  to  transport.  Still  another 
advantage  is  that  the  discrete  unknowns  are  defined  as  edge  or  face  fluxes,  which  are  finite 
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and  bounded  even  as  a  point  source  or  sink  is  approached,  so  that  the  MFEM  is  able  to 
handle  injection  and  pumping  wells  more  accurately  than  other  methods. 

Solution  of  discrete  equations.  A  major  obstacle  to  three-dimensional  applications  of  MFEM 
is  the  difficulty  of  solving  the  discrete  linear  equations  iteratively.  The  system  is  indefinite, 
and  the  usual  iterative  methods  are  unreliable.  The  quest  for  efficient,  reliable  solvers  has 
been  active  for  many  years.  A  successful  two-dimensional  algorithm  originated  by  Chavent 
et  al.  in  1984  uses  a  local  basis  for  the  subspace  of  divergence-free  functions  to  decompose  the 
problem  into  three  parts,  each  of  which  can  be  computed  efficiently.  Most  of  the  computation 
is  a  positive-definite  projection  into  this  subspace,  and  the  algorithm’s  efficiency  derives  in 
part  from  the  fact  that  the  subspace  dimension  is  substantially  less  than  that  of  the  full 
velocity  space. 

References  9  and  17  document  the  theory  behind  an  extension  of  this  solver  to  three 
dimensions.  Unlike  some  other  algorithms,  this  works  when  the  hydraulic  conductivity  is 
anisotropic,  resulting  in  cross-derivative  terms  that  create  additional  non-zero  coefficients 
in  the  linear  equations.  Numerical  results,  in  References  16,  17,  23,  and  31,  show  that  the 
convergence  rate  of  this  solver  is  independent  of  the  grid  size,  as  predicted  theoretically.  In 
a  pleasant  surprise,  the  convergence  rate  in  most  cases  is  also  independent  of  the  magnitude 
of  variations  in  hydraulic  conductivity. 

Distorted  grids.  Another  restriction  on  the  applicability  of  the  Raviart-Thomas  MFEM  de¬ 
scribed  above  is  that,  for  optimal  accuracy,  the  elements  must  be  rectangles,  triangles,  or 
parallelograms  (some  of  each  is  permissible).  Existing  finite  difference  grids  and  data  struc¬ 
tures  are  primarily  based  on  a  logically  rectangular  connectivity  pattern.  To  the  extent 
possible,  it  would  be  desirable  to  treat  distorted  logically  rectangular  grids,  of  quadrilater¬ 
als  in  two  dimensions  and  hexahedra  in  three,  designed  to  conform  to  the  irregularities  of 
subsurface  geology. 

Control-volume  mixed  finite  element  method.  With  this  in  mind,  this  project  designed  and 
implemented  in  two  (References  7,  11,  18,  30)  and  three  (References  13,  16,  20,  23,  31) 
dimensions  a  variant  of  MFEM,  called  the  control-volume  mixed  finite  element  method 
(CVMFEM),  that  preserves  optimal-order  accuracy  on  distorted  grids.  In  this  procedure, 
the  velocity  test  functions  are  altered  from  those  of  the  MFEM,  in  such  a  way  as  to  be 
distortions  of  piecewise  constants  on  control  volumes.  The  result  is  is  a  discrete  Darcy 
law  for  a  distorted  grid:  the  pressure  difference  between  two  adjoining  cells  is,  instead  of 
a  multiple  of  the  flux  across  their  common  edge  or  face,  a  linear  combination  of  the  fluxes 
across  all  of  their  edges  or  faces. 

The  reputation  of  MFEMs  is  that  they  provide  more  accurate  velocities  and  streamlines 
in  problems  with  heterogeneous  hydraulic  conductivity;  for  CVMFEM,  Reference  7  showed 
this  with  both  mathematical  theory  and  two-dimensional  numerical  computations  that  in¬ 
cluded  anisotropy  and  grid  distortions.  Other  methods  have  particular  difficulties  with  large 
discontinuous  jumps  in  conductivity,  which  are  of  obvious  importance  for  subsurface  flow,  but 
MFEMs  do  not.  CVMFEM  computed  fluxes  to  the  accuracy  of  the  square  of  the  grid  size, 
i.e.,  second-order  accuracy.  Subsequently,  similar  results  were  found  for  three-dimensional 
CVMFEM  in  References  13,  20,  and  31.  The  efficient  equation  solver  was  adapted  to  MFEM 
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and  CVMFEM  with  distorted  grids  (References  16,  23,  31).  The  essential  structure  of  the 
solver  is  the  same  as  for  rectangular  grids,  but  the  determination  of  the  equation  coeffi¬ 
cients  is  more  complicated.  A  CVMFEM-based  flow  code  is  currently  under  development 
in  collaboration  with  USGS,  and  it  may  eventually  become  a  successor  to  the  widely-used 
MODFLOW. 

Justification  of  model  equations 

An  important  issue  in  modeling  of  contaminant  transport  at  DOD  sites  is  justification  of 
the  manner  in  which  the  model  equations  represent  the  physical  processes.  Subsurface 
formations  are  heterogeneous  at  a  multitude  of  length  scales,  including  heterogeneities  at 
scales  finer  than  a  computational  grid.  The  effects  of  such  properties  must  be  represented 
indirectly  in  models,  through  such  concepts  as  dispersivity  in  solute  transport.  References 
3  and  10  report  on  a  formulation,  implemented  with  stochastic  finite  element  methods, 
designed  to  compute  effective  dispersivities  for  real  data.  Such  data  may  violate  the  idealized 
assumptions  that  underlie  analytical  approaches  to  this  problem.  Reference  2,  an  invited 
review  paper,  discusses  the  physical  assumptions  made  by  various  models  of  multiphase 
multicomponent  transport,  with  particular  attention  to  the  local  equilibrium  assumption 
and  the  implications  for  practical  modeling  when  it  is  valid  or  invalid. 

Other  research 

References  4  and  6  document  investigations  into  biased  interpretations  of  field  measurements. 
A  mathematical  model  consisting  of  a  diffusion  equation  with  a  discontinuous  piecewise- 
constant  coefficient,  together  with  analytical  and  numerical  solutions,  showed  that  the  usual 
calibrations  used  by  field  scientists  could  be  substantially  in  error  in  some  cases  of  practical 
interest.  Reference  27  studied  operator-splitting  techniques  for  reactive  transport,  delin¬ 
eating  types  of  cases  in  which  it  was  advantageous  or  disadvantageous  to  iterate,  based 
on  stability  analysis.  Reference  28  described  a  robust  incomplete-LU  (ILU)  preconditioner 
for  conjugate- gradient  solution  of  linear  equations  arising  from  flow  problems  with  varying 
anisotropy  and  heterogeneity.  The  idea  was  to  minimize  the  directional  sensitivity  of  the 
speed  of  convergence  by  performing  the  ILU  decomposition  in  a  zigzag  fashion  through  the 
grid. 
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